#!/usr/bin/env julia

"""
    mock_data_generator.jl

Comprehensive mock data generation utilities for testing all observation processing modules.
This module provides realistic synthetic observation data for comprehensive testing of
data assimilation algorithms while ensuring deterministic and reproducible test results.

# Mock Data Categories

## Observation Data
- Conventional observations (surface, radiosonde, aircraft, marine)
- Satellite observations (radiance, winds, retrievals)
- Specialized observations (GPS RO, radar, lightning, aerosol)

## Atmospheric Fields
- Background state variables (temperature, pressure, humidity, winds)
- Error covariance matrices (B matrix components)
- Model grids and coordinate systems

## Quality Control Scenarios
- Good quality observations
- Poor quality observations with various error types
- Outliers and gross errors
- Missing data patterns

# Design Principles

- **Deterministic**: All random generation uses fixed seeds
- **Realistic**: Physically consistent values and error characteristics
- **Comprehensive**: Coverage of all observation types and scenarios
- **Scalable**: Efficient generation for large test datasets
- **Configurable**: Easy to modify for different test scenarios
"""

module MockDataGenerator

using Random
using Dates
using Statistics
using LinearAlgebra
using StaticArrays

# Use relative imports to avoid circular dependencies
# Import observation types from the main module when available
using Test  # For @test_nowarn in generators

export MockConfig, generate_mock_observations, generate_mock_background
export generate_conventional_observations, generate_satellite_observations
export generate_specialized_observations, generate_atmospheric_field
export create_test_grid, generate_error_covariances
export MockQualityScenario, apply_quality_scenario

"""Configuration for mock data generation"""
struct MockConfig{T<:AbstractFloat}
    # Spatial domain
    lat_range::Tuple{T, T}
    lon_range::Tuple{T, T}
    pressure_levels::Vector{T}
    
    # Temporal domain
    time_start::DateTime
    time_end::DateTime
    time_interval::Minute
    
    # Observation counts
    n_surface::Int
    n_radiosonde::Int
    n_aircraft::Int
    n_satellite_ir::Int
    n_satellite_mw::Int
    n_gps_ro::Int
    n_radar::Int
    
    # Error characteristics
    typical_temp_error::T
    typical_pressure_error::T
    typical_humidity_error::T
    typical_wind_error::T
    typical_radiance_error::T
    
    # Quality control parameters
    good_data_fraction::T
    outlier_fraction::T
    missing_data_fraction::T
    
    # Random seed for reproducibility
    random_seed::Int
    
    function MockConfig{T}(;
        lat_range::Tuple{T, T} = (-90.0, 90.0),
        lon_range::Tuple{T, T} = (-180.0, 180.0),
        pressure_levels::Vector{T} = T[100000, 85000, 70000, 50000, 30000, 20000, 10000],
        time_start::DateTime = DateTime(2024, 1, 1, 0, 0, 0),
        time_end::DateTime = DateTime(2024, 1, 1, 6, 0, 0),
        time_interval::Minute = Minute(60),
        n_surface::Int = 1000,
        n_radiosonde::Int = 200,
        n_aircraft::Int = 500,
        n_satellite_ir::Int = 10000,
        n_satellite_mw::Int = 5000,
        n_gps_ro::Int = 100,
        n_radar::Int = 50,
        typical_temp_error::T = 1.0,
        typical_pressure_error::T = 100.0,
        typical_humidity_error::T = 0.1,
        typical_wind_error::T = 2.0,
        typical_radiance_error::T = 0.5,
        good_data_fraction::T = 0.85,
        outlier_fraction::T = 0.05,
        missing_data_fraction::T = 0.10,
        random_seed::Int = 42
    ) where T
        new(lat_range, lon_range, pressure_levels, time_start, time_end, time_interval,
            n_surface, n_radiosonde, n_aircraft, n_satellite_ir, n_satellite_mw,
            n_gps_ro, n_radar, typical_temp_error, typical_pressure_error,
            typical_humidity_error, typical_wind_error, typical_radiance_error,
            good_data_fraction, outlier_fraction, missing_data_fraction, random_seed)
    end
end

# Convenience constructor for Float64
MockConfig(; kwargs...) = MockConfig{Float64}(; kwargs...)

"""Quality control scenario types for testing"""
@enum MockQualityScenario begin
    GOOD_QUALITY
    POOR_QUALITY
    GROSS_ERRORS
    MISSING_DATA
    MIXED_QUALITY
end

"""
    generate_mock_observations(config::MockConfig{T}) -> Dict{String, Vector}

Generate comprehensive mock observation dataset for testing.
"""
function generate_mock_observations(config::MockConfig{T}) where T
    Random.seed!(config.random_seed)
    
    observations = Dict{String, Vector}()
    
    # Generate conventional observations
    observations["surface"] = generate_surface_observations(config)
    observations["radiosonde"] = generate_radiosonde_observations(config)
    observations["aircraft"] = generate_aircraft_observations(config)
    observations["marine"] = generate_marine_observations(config)
    
    # Generate satellite observations
    observations["satellite_ir"] = generate_satellite_ir_observations(config)
    observations["satellite_mw"] = generate_satellite_mw_observations(config)
    observations["satellite_winds"] = generate_satellite_wind_observations(config)
    
    # Generate specialized observations
    observations["gps_ro"] = generate_gps_ro_observations(config)
    observations["radar"] = generate_radar_observations(config)
    observations["lightning"] = generate_lightning_observations(config)
    
    return observations
end

"""Generate realistic surface observations"""
function generate_surface_observations(config::MockConfig{T}) where T
    observations = []
    
    for i in 1:config.n_surface
        # Realistic geographic distribution (more observations over land)
        lat = rand() * (config.lat_range[2] - config.lat_range[1]) + config.lat_range[1]
        lon = rand() * (config.lon_range[2] - config.lon_range[1]) + config.lon_range[1]
        
        # Realistic elevation distribution
        elevation = max(0.0, randn() * 500.0 + 200.0)  # Sea level to mountains
        
        # Generate realistic meteorological values
        # Temperature varies with latitude and elevation
        base_temp = 288.15 - 0.0065 * elevation  # Standard atmosphere lapse rate
        seasonal_temp = 15.0 * sin(2π * (lat + 90) / 180)  # Seasonal variation
        temp = base_temp + seasonal_temp + randn() * 5.0
        
        # Surface pressure varies with elevation
        surface_pressure = 101325.0 * exp(-0.0001225 * elevation)
        surface_pressure += randn() * config.typical_pressure_error
        
        # Humidity varies with temperature and latitude
        relative_humidity = min(100.0, max(0.0, 60.0 + 30.0 * sin(π * lat / 180) + randn() * 20.0))
        
        # Wind components
        wind_u = randn() * 5.0 + 2.0 * sin(π * lat / 90)  # Westerlies
        wind_v = randn() * 3.0
        
        # Observation time
        time_range = Dates.value(config.time_end - config.time_start) / 1000  # milliseconds
        obs_time = config.time_start + Millisecond(rand() * time_range)
        
        # Create observation dictionary (representing the actual observation type)
        obs = Dict{String, Any}(
            "obs_type" => "surface",
            "variables" => Dict(
                "temperature" => temp,
                "surface_pressure" => surface_pressure,
                "relative_humidity" => relative_humidity,
                "wind_u" => wind_u,
                "wind_v" => wind_v
            ),
            "errors" => Dict(
                "temperature" => config.typical_temp_error * (0.8 + 0.4 * rand()),
                "surface_pressure" => config.typical_pressure_error * (0.8 + 0.4 * rand()),
                "relative_humidity" => 5.0 * (0.8 + 0.4 * rand()),
                "wind_u" => config.typical_wind_error * (0.8 + 0.4 * rand()),
                "wind_v" => config.typical_wind_error * (0.8 + 0.4 * rand())
            ),
            "location" => Dict(
                "latitude" => lat,
                "longitude" => lon,
                "elevation" => elevation
            ),
            "time" => obs_time,
            "qc_flags" => Dict{String, Int}(),
            "metadata" => Dict(
                "station_id" => "MOCK_$(lpad(i, 6, '0'))",
                "station_type" => "land"
            )
        )
        
        push!(observations, obs)
    end
    
    return observations
end

"""Generate realistic radiosonde observations"""
function generate_radiosonde_observations(config::MockConfig{T}) where T
    observations = []
    
    for i in 1:config.n_radiosonde
        # Radiosonde locations (typically at weather stations)
        lat = rand() * (config.lat_range[2] - config.lat_range[1]) + config.lat_range[1]
        lon = rand() * (config.lon_range[2] - config.lon_range[1]) + config.lon_range[1]
        elevation = max(0.0, randn() * 300.0 + 100.0)
        
        # Standard radiosonde launch times (00Z, 12Z typically)
        launch_hour = rand([0, 6, 12, 18])
        obs_time = DateTime(2024, 1, 1, launch_hour, 0, 0)
        
        # Generate vertical profile
        profile = []
        for p_level in config.pressure_levels
            if p_level >= 50000  # Only generate data for lower atmosphere in this example
                # Temperature profile with realistic lapse rate
                altitude = -7000 * log(p_level / 101325.0)  # Approximate altitude
                temp = 288.15 - 0.0065 * altitude + randn() * 2.0
                
                # Humidity decreases with altitude
                rel_humidity = min(100.0, max(0.0, 80.0 * exp(-altitude / 8000) + randn() * 15.0))
                
                # Wind increases with altitude (jet stream effect)
                wind_u = 10.0 + 20.0 * exp(-(p_level - 25000)^2 / (2 * 15000^2)) + randn() * 3.0
                wind_v = randn() * 5.0
                
                level_data = Dict(
                    "pressure" => p_level,
                    "temperature" => temp,
                    "relative_humidity" => rel_humidity,
                    "wind_u" => wind_u,
                    "wind_v" => wind_v,
                    "altitude" => altitude
                )
                
                push!(profile, level_data)
            end
        end
        
        # Create radiosonde observation
        obs = Dict{String, Any}(
            "obs_type" => "radiosonde",
            "profile" => profile,
            "location" => Dict(
                "latitude" => lat,
                "longitude" => lon,
                "elevation" => elevation
            ),
            "time" => obs_time,
            "qc_flags" => Dict{String, Int}(),
            "metadata" => Dict(
                "station_id" => "RS_$(lpad(i, 4, '0'))",
                "balloon_type" => "rubber",
                "max_altitude" => maximum([level["altitude"] for level in profile])
            )
        )
        
        push!(observations, obs)
    end
    
    return observations
end

"""Generate realistic aircraft observations"""
function generate_aircraft_observations(config::MockConfig{T}) where T
    observations = []
    
    for i in 1:config.n_aircraft
        # Aircraft typically fly at commercial flight levels
        flight_level = rand([250, 300, 350, 390, 410, 430]) * 100  # hPa equivalent
        altitude = -7000 * log(flight_level / 101325.0)
        
        # Flight paths typically follow great circles
        lat = rand() * (config.lat_range[2] - config.lat_range[1]) + config.lat_range[1]
        lon = rand() * (config.lon_range[2] - config.lon_range[1]) + config.lon_range[1]
        
        # Atmospheric conditions at flight level
        temp = 218.15 + randn() * 5.0  # Typical tropopause temperature
        wind_u = 30.0 + randn() * 15.0  # Strong winds at flight level
        wind_v = randn() * 10.0
        
        obs_time = config.time_start + Millisecond(rand() * 
                   Dates.value(config.time_end - config.time_start))
        
        obs = Dict{String, Any}(
            "obs_type" => "aircraft",
            "variables" => Dict(
                "temperature" => temp,
                "wind_u" => wind_u,
                "wind_v" => wind_v
            ),
            "errors" => Dict(
                "temperature" => config.typical_temp_error * 1.2,  # Slightly higher error
                "wind_u" => config.typical_wind_error * 1.5,
                "wind_v" => config.typical_wind_error * 1.5
            ),
            "location" => Dict(
                "latitude" => lat,
                "longitude" => lon,
                "altitude" => altitude,
                "pressure" => flight_level
            ),
            "time" => obs_time,
            "qc_flags" => Dict{String, Int}(),
            "metadata" => Dict(
                "aircraft_id" => "AC_$(lpad(i, 5, '0'))",
                "flight_phase" => rand(["cruise", "ascent", "descent"])
            )
        )
        
        push!(observations, obs)
    end
    
    return observations
end

"""Generate marine observations (ships, buoys)"""
function generate_marine_observations(config::MockConfig{T}) where T
    observations = []
    n_marine = config.n_surface ÷ 10  # About 10% of surface obs are marine
    
    for i in 1:n_marine
        # Marine observations over oceans/seas
        lat = rand() * (config.lat_range[2] - config.lat_range[1]) + config.lat_range[1]
        lon = rand() * (config.lon_range[2] - config.lon_range[1]) + config.lon_range[1]
        
        # Sea surface temperature
        sst = 273.15 + 20.0 + 10.0 * sin(π * lat / 90) + randn() * 3.0
        
        # Marine atmospheric conditions
        temp = sst + randn() * 2.0  # Air temp close to SST
        surface_pressure = 101325.0 + randn() * config.typical_pressure_error
        rel_humidity = min(100.0, max(50.0, 85.0 + randn() * 10.0))  # Higher humidity over ocean
        
        # Ocean winds
        wind_u = randn() * 8.0
        wind_v = randn() * 8.0
        
        obs_time = config.time_start + Millisecond(rand() * 
                   Dates.value(config.time_end - config.time_start))
        
        obs = Dict{String, Any}(
            "obs_type" => "marine",
            "variables" => Dict(
                "temperature" => temp,
                "surface_pressure" => surface_pressure,
                "relative_humidity" => rel_humidity,
                "wind_u" => wind_u,
                "wind_v" => wind_v,
                "sea_surface_temperature" => sst
            ),
            "errors" => Dict(
                "temperature" => config.typical_temp_error * 1.1,
                "surface_pressure" => config.typical_pressure_error * 0.9,
                "relative_humidity" => 8.0,
                "wind_u" => config.typical_wind_error * 1.3,
                "wind_v" => config.typical_wind_error * 1.3,
                "sea_surface_temperature" => 0.5
            ),
            "location" => Dict(
                "latitude" => lat,
                "longitude" => lon,
                "elevation" => 0.0  # Sea level
            ),
            "time" => obs_time,
            "qc_flags" => Dict{String, Int}(),
            "metadata" => Dict(
                "platform_id" => "MARINE_$(lpad(i, 4, '0'))",
                "platform_type" => rand(["ship", "buoy", "oil_platform"])
            )
        )
        
        push!(observations, obs)
    end
    
    return observations
end

"""Generate satellite infrared radiance observations"""
function generate_satellite_ir_observations(config::MockConfig{T}) where T
    observations = []
    
    # Typical IR channels for GOES, SEVIRI, etc.
    ir_channels = [6.2, 6.9, 7.3, 8.5, 9.6, 10.4, 11.2, 12.0, 13.3]  # μm
    
    for i in 1:config.n_satellite_ir
        # Satellite viewing geometry
        lat = rand() * (config.lat_range[2] - config.lat_range[1]) + config.lat_range[1]
        lon = rand() * (config.lon_range[2] - config.lon_range[1]) + config.lon_range[1]
        
        satellite_zenith = rand() * 60.0  # degrees
        satellite_azimuth = rand() * 360.0
        
        # Select random channel
        channel = rand(ir_channels)
        
        # Brightness temperature varies by channel and atmospheric conditions
        if channel < 7.0  # Water vapor channels
            bt = 200.0 + 50.0 * rand()
        elseif channel < 10.0  # Mid-level channels
            bt = 220.0 + 40.0 * rand()
        else  # Window channels
            bt = 250.0 + 30.0 * rand()
        end
        
        bt += randn() * config.typical_radiance_error
        
        obs_time = config.time_start + Millisecond(rand() * 
                   Dates.value(config.time_end - config.time_start))
        
        obs = Dict{String, Any}(
            "obs_type" => "satellite_radiance",
            "sensor_type" => "infrared",
            "variables" => Dict(
                "brightness_temperature" => bt
            ),
            "errors" => Dict(
                "brightness_temperature" => config.typical_radiance_error * (0.8 + 0.4 * rand())
            ),
            "location" => Dict(
                "latitude" => lat,
                "longitude" => lon,
                "satellite_zenith_angle" => satellite_zenith,
                "satellite_azimuth_angle" => satellite_azimuth
            ),
            "time" => obs_time,
            "channel_info" => Dict(
                "wavelength" => channel,
                "channel_number" => findfirst(x -> x == channel, ir_channels)
            ),
            "qc_flags" => Dict{String, Int}(),
            "metadata" => Dict(
                "satellite_id" => rand(["GOES-16", "GOES-17", "MSG-4", "Himawari-8"]),
                "sensor_id" => "IR_imager"
            )
        )
        
        push!(observations, obs)
    end
    
    return observations
end

"""Generate satellite microwave radiance observations"""
function generate_satellite_mw_observations(config::MockConfig{T}) where T
    observations = []
    
    # AMSU-A like channels (K)
    amsu_frequencies = [23.8, 31.4, 50.3, 52.8, 53.6, 54.4, 54.9, 55.5, 57.3, 58.0]  # GHz
    
    for i in 1:config.n_satellite_mw
        lat = rand() * (config.lat_range[2] - config.lat_range[1]) + config.lat_range[1]
        lon = rand() * (config.lon_range[2] - config.lon_range[1]) + config.lon_range[1]
        
        satellite_zenith = rand() * 45.0  # Microwave sensors have narrower swath
        scan_position = rand(1:30)  # Scan position across swath
        
        # Select random channel
        frequency = rand(amsu_frequencies)
        
        # Brightness temperature varies by frequency
        if frequency < 30.0  # Window channels
            bt = 200.0 + 80.0 * rand()
        elseif frequency < 55.0  # Temperature sounding channels
            bt = 150.0 + 100.0 * rand()
        else  # Upper atmosphere channels
            bt = 100.0 + 100.0 * rand()
        end
        
        bt += randn() * config.typical_radiance_error * 1.5  # Higher MW error
        
        obs_time = config.time_start + Millisecond(rand() * 
                   Dates.value(config.time_end - config.time_start))
        
        obs = Dict{String, Any}(
            "obs_type" => "satellite_radiance",
            "sensor_type" => "microwave",
            "variables" => Dict(
                "brightness_temperature" => bt
            ),
            "errors" => Dict(
                "brightness_temperature" => config.typical_radiance_error * 1.5 * (0.8 + 0.4 * rand())
            ),
            "location" => Dict(
                "latitude" => lat,
                "longitude" => lon,
                "satellite_zenith_angle" => satellite_zenith,
                "scan_position" => scan_position
            ),
            "time" => obs_time,
            "channel_info" => Dict(
                "frequency" => frequency,
                "channel_number" => findfirst(x -> x == frequency, amsu_frequencies)
            ),
            "qc_flags" => Dict{String, Int}(),
            "metadata" => Dict(
                "satellite_id" => rand(["NOAA-20", "NOAA-21", "Metop-C"]),
                "sensor_id" => "AMSU-A"
            )
        )
        
        push!(observations, obs)
    end
    
    return observations
end

"""Generate satellite atmospheric motion vector (wind) observations"""
function generate_satellite_wind_observations(config::MockConfig{T}) where T
    observations = []
    n_winds = config.n_satellite_ir ÷ 5  # Fewer wind observations than radiances
    
    for i in 1:n_winds
        lat = rand() * (config.lat_range[2] - config.lat_range[1]) + config.lat_range[1]
        lon = rand() * (config.lon_range[2] - config.lon_range[1]) + config.lon_range[1]
        
        # Wind assignment pressure level
        pressure_level = rand([200, 400, 700, 850, 925]) * 100.0  # Pa
        
        # Realistic wind speeds vary by level
        if pressure_level < 30000  # Upper level
            wind_speed = 20.0 + 30.0 * rand()
        elseif pressure_level < 70000  # Mid level
            wind_speed = 10.0 + 20.0 * rand()
        else  # Low level
            wind_speed = 5.0 + 15.0 * rand()
        end
        
        wind_direction = rand() * 360.0  # degrees
        
        # Convert to u, v components
        wind_u = wind_speed * sin(deg2rad(wind_direction))
        wind_v = wind_speed * cos(deg2rad(wind_direction))
        
        obs_time = config.time_start + Millisecond(rand() * 
                   Dates.value(config.time_end - config.time_start))
        
        obs = Dict{String, Any}(
            "obs_type" => "satellite_winds",
            "variables" => Dict(
                "wind_u" => wind_u,
                "wind_v" => wind_v,
                "wind_speed" => wind_speed,
                "wind_direction" => wind_direction
            ),
            "errors" => Dict(
                "wind_u" => config.typical_wind_error * 1.8,  # Higher error for satellite winds
                "wind_v" => config.typical_wind_error * 1.8
            ),
            "location" => Dict(
                "latitude" => lat,
                "longitude" => lon,
                "pressure" => pressure_level
            ),
            "time" => obs_time,
            "qc_flags" => Dict{String, Int}(),
            "metadata" => Dict(
                "satellite_id" => rand(["GOES-16", "GOES-17", "Meteosat-11"]),
                "derivation_method" => "IR_tracer",
                "quality_indicator" => rand(["good", "fair", "poor"])
            )
        )
        
        push!(observations, obs)
    end
    
    return observations
end

"""Generate GPS Radio Occultation observations"""
function generate_gps_ro_observations(config::MockConfig{T}) where T
    observations = []
    
    for i in 1:config.n_gps_ro
        lat = rand() * (config.lat_range[2] - config.lat_range[1]) + config.lat_range[1]
        lon = rand() * (config.lon_range[2] - config.lon_range[1]) + config.lon_range[1]
        
        # GPS RO provides vertical profiles
        profile = []
        altitudes = collect(5000:1000:30000)  # 5-30 km altitude range
        
        for alt in altitudes
            # Refractivity decreases exponentially with altitude
            refractivity = 300.0 * exp(-alt / 7000.0) * (0.8 + 0.4 * rand())
            
            # Bending angle
            impact_parameter = 6371000.0 + alt  # Earth radius + altitude
            bending_angle = refractivity * 1e-6 + randn() * 5e-6
            
            level_data = Dict(
                "altitude" => alt,
                "impact_parameter" => impact_parameter,
                "refractivity" => refractivity,
                "bending_angle" => bending_angle
            )
            
            push!(profile, level_data)
        end
        
        obs_time = config.time_start + Millisecond(rand() * 
                   Dates.value(config.time_end - config.time_start))
        
        obs = Dict{String, Any}(
            "obs_type" => "gps_ro",
            "profile" => profile,
            "location" => Dict(
                "latitude" => lat,
                "longitude" => lon
            ),
            "time" => obs_time,
            "qc_flags" => Dict{String, Int}(),
            "metadata" => Dict(
                "leo_satellite" => rand(["COSMIC-2", "Spire", "MetOp"]),
                "gps_satellite" => "GPS_$(rand(1:32))",
                "occultation_id" => "RO_$(lpad(i, 6, '0'))"
            )
        )
        
        push!(observations, obs)
    end
    
    return observations
end

"""Generate radar observations (precipitation, winds)"""
function generate_radar_observations(config::MockConfig{T}) where T
    observations = []
    
    for i in 1:config.n_radar
        # Radar locations (typically fixed sites)
        lat = rand() * (config.lat_range[2] - config.lat_range[1]) + config.lat_range[1]
        lon = rand() * (config.lon_range[2] - config.lon_range[1]) + config.lon_range[1]
        
        # Radar scan geometry
        range_km = 50.0 + 200.0 * rand()  # 50-250 km range
        azimuth = rand() * 360.0
        elevation_angle = 0.5 + 4.5 * rand()  # 0.5-5.0 degrees
        
        # Radar measurements
        reflectivity = -20.0 + 60.0 * rand()  # dBZ
        radial_velocity = -30.0 + 60.0 * rand()  # m/s
        
        obs_time = config.time_start + Millisecond(rand() * 
                   Dates.value(config.time_end - config.time_start))
        
        obs = Dict{String, Any}(
            "obs_type" => "radar",
            "variables" => Dict(
                "reflectivity" => reflectivity,
                "radial_velocity" => radial_velocity
            ),
            "errors" => Dict(
                "reflectivity" => 2.0,  # dBZ
                "radial_velocity" => 1.0  # m/s
            ),
            "location" => Dict(
                "radar_latitude" => lat,
                "radar_longitude" => lon,
                "target_range" => range_km * 1000.0,  # meters
                "azimuth" => azimuth,
                "elevation" => elevation_angle
            ),
            "time" => obs_time,
            "qc_flags" => Dict{String, Int}(),
            "metadata" => Dict(
                "radar_id" => "RADAR_$(lpad(i, 3, '0'))",
                "radar_type" => "WSR-88D",
                "scan_mode" => "precipitation"
            )
        )
        
        push!(observations, obs)
    end
    
    return observations
end

"""Generate lightning observations"""
function generate_lightning_observations(config::MockConfig{T}) where T
    observations = []
    n_lightning = config.n_radar ÷ 2  # Fewer lightning observations
    
    for i in 1:n_lightning
        lat = rand() * (config.lat_range[2] - config.lat_range[1]) + config.lat_range[1]
        lon = rand() * (config.lon_range[2] - config.lon_range[1]) + config.lon_range[1]
        
        # Lightning characteristics
        peak_current = 5000.0 + 45000.0 * rand()  # Amperes
        multiplicity = rand(1:8)  # Number of strokes
        
        obs_time = config.time_start + Millisecond(rand() * 
                   Dates.value(config.time_end - config.time_start))
        
        obs = Dict{String, Any}(
            "obs_type" => "lightning",
            "variables" => Dict(
                "peak_current" => peak_current,
                "multiplicity" => multiplicity
            ),
            "errors" => Dict(
                "peak_current" => peak_current * 0.2,  # 20% error
                "location" => 2000.0  # 2 km location error
            ),
            "location" => Dict(
                "latitude" => lat,
                "longitude" => lon
            ),
            "time" => obs_time,
            "qc_flags" => Dict{String, Int}(),
            "metadata" => Dict(
                "network" => rand(["NLDN", "GLD360", "WWLLN"]),
                "stroke_type" => rand(["cloud-to-ground", "intracloud"])
            )
        )
        
        push!(observations, obs)
    end
    
    return observations
end

"""
    apply_quality_scenario!(observations::Vector, scenario::MockQualityScenario)

Apply quality control scenarios to mock observations for testing QC algorithms.
"""
function apply_quality_scenario!(observations::Vector, scenario::MockQualityScenario)
    if scenario == GOOD_QUALITY
        # All observations remain good quality (default)
        for obs in observations
            obs["qc_flags"]["gross_error"] = 0
            obs["qc_flags"]["background_check"] = 0
            obs["qc_flags"]["buddy_check"] = 0
        end
        
    elseif scenario == POOR_QUALITY
        # Add random noise and increase error estimates
        for obs in observations
            if haskey(obs, "variables")
                for (var, value) in obs["variables"]
                    if typeof(value) <: Real
                        noise_factor = 1.0 + randn() * 0.5  # Up to 50% noise
                        obs["variables"][var] = value * noise_factor
                        if haskey(obs["errors"], var)
                            obs["errors"][var] *= 2.0  # Double error estimates
                        end
                    end
                end
            end
            obs["qc_flags"]["gross_error"] = 0
            obs["qc_flags"]["background_check"] = 1  # Flag as questionable
        end
        
    elseif scenario == GROSS_ERRORS
        # Add gross errors to subset of observations
        n_gross = max(1, length(observations) ÷ 10)  # 10% gross errors
        gross_indices = rand(1:length(observations), n_gross)
        
        for (i, obs) in enumerate(observations)
            if i in gross_indices
                # Add gross errors
                if haskey(obs, "variables")
                    for (var, value) in obs["variables"]
                        if typeof(value) <: Real
                            error_magnitude = 5.0 + 10.0 * rand()  # Large error
                            error_sign = rand([-1, 1])
                            obs["variables"][var] = value + error_sign * error_magnitude * 
                                                   get(obs["errors"], var, 1.0)
                        end
                    end
                end
                obs["qc_flags"]["gross_error"] = 1
            else
                obs["qc_flags"]["gross_error"] = 0
            end
        end
        
    elseif scenario == MISSING_DATA
        # Remove random subset of observations or variables
        n_missing = max(1, length(observations) ÷ 5)  # 20% missing
        missing_indices = rand(1:length(observations), n_missing)
        
        for (i, obs) in enumerate(observations)
            if i in missing_indices
                obs["qc_flags"]["missing_data"] = 1
                # Randomly remove some variables
                if haskey(obs, "variables") && length(obs["variables"]) > 1
                    vars_to_keep = rand(1:length(obs["variables"])-1)
                    var_keys = collect(keys(obs["variables"]))
                    keep_keys = var_keys[1:vars_to_keep]
                    new_vars = Dict()
                    new_errors = Dict()
                    for key in keep_keys
                        new_vars[key] = obs["variables"][key]
                        if haskey(obs["errors"], key)
                            new_errors[key] = obs["errors"][key]
                        end
                    end
                    obs["variables"] = new_vars
                    obs["errors"] = new_errors
                end
            else
                obs["qc_flags"]["missing_data"] = 0
            end
        end
        
    elseif scenario == MIXED_QUALITY
        # Apply mixed quality control scenarios
        n_obs = length(observations)
        n_good = Int(round(0.7 * n_obs))
        n_poor = Int(round(0.2 * n_obs))
        n_gross = Int(round(0.1 * n_obs))
        
        # Shuffle indices and assign quality types
        indices = shuffle(1:n_obs)
        good_indices = Set(indices[1:n_good])
        poor_indices = Set(indices[n_good+1:n_good+n_poor])
        gross_indices = Set(indices[n_good+n_poor+1:end])
        
        for (i, obs) in enumerate(observations)
            if i in good_indices
                obs["qc_flags"]["gross_error"] = 0
                obs["qc_flags"]["background_check"] = 0
                obs["qc_flags"]["buddy_check"] = 0
            elseif i in poor_indices
                # Add moderate noise
                if haskey(obs, "variables")
                    for (var, value) in obs["variables"]
                        if typeof(value) <: Real
                            noise = randn() * get(obs["errors"], var, 1.0) * 0.8
                            obs["variables"][var] = value + noise
                        end
                    end
                end
                obs["qc_flags"]["gross_error"] = 0
                obs["qc_flags"]["background_check"] = 1
            else  # gross error
                # Add large errors
                if haskey(obs, "variables")
                    for (var, value) in obs["variables"]
                        if typeof(value) <: Real
                            error_magnitude = 3.0 + 7.0 * rand()
                            error_sign = rand([-1, 1])
                            obs["variables"][var] = value + error_sign * error_magnitude * 
                                                   get(obs["errors"], var, 1.0)
                        end
                    end
                end
                obs["qc_flags"]["gross_error"] = 1
                obs["qc_flags"]["background_check"] = 1
            end
        end
    end
    
    return observations
end

"""
    create_test_grid(nx::Int, ny::Int, nz::Int) -> Dict

Create a test atmospheric grid for forward operator testing.
"""
function create_test_grid(nx::Int=32, ny::Int=24, nz::Int=16)
    # Create regular lat-lon grid
    lats = range(-90.0, 90.0, length=ny)
    lons = range(-180.0, 180.0, length=nx)
    
    # Pressure levels (Pa)
    pressure_levels = [
        100000, 92500, 85000, 77500, 70000, 60000, 50000, 40000,
        30000, 25000, 20000, 15000, 10000, 7000, 5000, 3000
    ][1:nz]
    
    return Dict(
        "latitude" => collect(lats),
        "longitude" => collect(lons),
        "pressure" => pressure_levels,
        "dimensions" => (nx, ny, nz)
    )
end

"""
    generate_atmospheric_field(grid::Dict, variable::String; smooth::Bool=true) -> Array{Float64, 3}

Generate realistic atmospheric field on the given grid.
"""
function generate_atmospheric_field(grid::Dict, variable::String; smooth::Bool=true)
    nx, ny, nz = grid["dimensions"]
    lats = grid["latitude"]
    lons = grid["longitude"]
    pressures = grid["pressure"]
    
    field = zeros(Float64, nx, ny, nz)
    
    for k in 1:nz, j in 1:ny, i in 1:nx
        lat = lats[j]
        lon = lons[i]
        p = pressures[k]
        
        if variable == "temperature"
            # Realistic temperature field with meridional gradient and vertical structure
            temp_surface = 288.15 + 20.0 * sin(π * lat / 180)  # Meridional temperature gradient
            altitude = -7000 * log(p / 101325.0)  # Approximate altitude
            field[i, j, k] = temp_surface - 0.0065 * altitude
            
            # Add some realistic variations
            field[i, j, k] += 5.0 * sin(2π * lon / 360) * exp(-altitude / 10000)  # Wave pattern
            
        elseif variable == "pressure"
            field[i, j, k] = p
            
        elseif variable == "humidity"
            # Humidity decreases with altitude and varies with temperature
            altitude = -7000 * log(p / 101325.0)
            humidity_surface = 0.015 * exp(-abs(lat) / 30)  # More humid in tropics
            field[i, j, k] = humidity_surface * exp(-altitude / 8000)
            
        elseif variable == "wind_u"
            # Westerly jet stream
            altitude = -7000 * log(p / 101325.0)
            jet_strength = 40.0 * exp(-(altitude - 10000)^2 / (2 * 3000^2))  # Jet at ~10km
            field[i, j, k] = jet_strength * exp(-(lat - 45)^2 / (2 * 15^2))  # Mid-latitude jet
            
        elseif variable == "wind_v"
            # Meridional wind component (much weaker)
            field[i, j, k] = 5.0 * sin(2π * lon / 360) * sin(π * lat / 180)
            
        else
            # Default: random field
            field[i, j, k] = randn()
        end
    end
    
    # Apply smoothing if requested
    if smooth && nx > 3 && ny > 3
        # Simple 3x3 smoothing kernel
        kernel = ones(3, 3) / 9.0
        for k in 1:nz
            # Apply 2D convolution (simplified)
            smoothed_slice = copy(field[:, :, k])
            for j in 2:ny-1, i in 2:nx-1
                smoothed_slice[i, j] = sum(field[i-1:i+1, j-1:j+1, k] .* kernel)
            end
            field[:, :, k] = smoothed_slice
        end
    end
    
    return field
end

"""
    generate_error_covariances(config::MockConfig{T}) -> Dict

Generate mock error covariance matrices for testing.
"""
function generate_error_covariances(config::MockConfig{T}) where T
    covariances = Dict{String, Matrix{T}}()
    
    # Observation error covariances (diagonal for simplicity)
    obs_types = ["temperature", "pressure", "humidity", "wind_u", "wind_v", "radiance"]
    typical_errors = [config.typical_temp_error, config.typical_pressure_error,
                     config.typical_humidity_error, config.typical_wind_error,
                     config.typical_wind_error, config.typical_radiance_error]
    
    for (obs_type, typical_error) in zip(obs_types, typical_errors)
        n_obs = 100  # Example size
        R = Matrix{T}(I, n_obs, n_obs) * typical_error^2
        
        # Add some off-diagonal correlations for realism
        for i in 1:n_obs-1
            correlation = 0.1 * exp(-abs(i - (i+1)) / 5.0)  # Short-range correlation
            R[i, i+1] = R[i+1, i] = correlation * typical_error^2
        end
        
        covariances["R_$(obs_type)"] = R
    end
    
    # Background error covariance (simplified)
    n_state = 1000  # Example state vector size
    B = Matrix{T}(I, n_state, n_state)
    
    # Add realistic correlations
    for i in 1:n_state
        for j in 1:n_state
            if i != j
                distance = abs(i - j)
                correlation = exp(-distance / 50.0)  # Exponential correlation
                B[i, j] = correlation * sqrt(B[i, i] * B[j, j])
            end
        end
    end
    
    covariances["B"] = B
    
    return covariances
end

end  # module MockDataGenerator